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ABSTRACT 

We present a time-dependent model of the interaction between a neutron star mag- 
netosphere and a thin (Shakura-Sunyaev) accretion disk, where the extent of the 
magnetosphere is determined by balancing outward diffusion and inward advection of 
the stellar magnetic field at the inner edge of the disc. The nature of the equilibria 
available to the system is governed by the magnetic Prandtl number Pm and the ratio 
£ of the corotation radius to the Alfven radius. For £ > Pm ' 3 , the system can occupy 
one of two stable states where the torques are of opposite signs. If the star is spinning 
up initially, in the absence of extraneous perturbations, £ decreases until the spin-up 
equilibrium vanishes, the star subsequently spins down, and the torque asymptotes to 
zero. Vortex-in-cell simulations of the Kelvin-Helmholtz instability suggest that the 
transport speed across the mixing layer between the disc and magnetosphere is less 
than the shear speed when the layer is thin, unlike in previous models. 

Key words: accretion, accretion discs - instabilities - MHD - pulsars: general - 
stars: neutron - X-rays: binaries 
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A complete time dependent theory of mass and angu- 
lar momentum transfer onto a neutron star accreting 
via a thin (Shakura-Sunyaev) disc needs to explain the 
stochastic fluctuations observed in the spin frequencies and 
X-ray luminosities of suc h systems feildsten et alJ Il997t 
iBavkal fc Oeeelmanl ll993Th as well as regular phenomena 
like the neg ative (spin-down) torque on accreting millisec - 
ond pulsars (iGallowav et al . 2002:; IChakrabartv et all2003t) . 
quasi-periodic os cillations (QPOs) in low-mass X-ray bi- 
naries (LMXBs) l|van der Klis 200 1|), and torque reversals 
jNelson et allll99itl . Torque reversals, in which some X- 
ray pulsars alternate between sustained episodes of spin- 
up/down, present a particular challenge to theories of accre- 
tion because the transitions between spin-up/down episodes 
are sudden (~ 10% of the duration of the episodes them- 
selves). 

In the standard iGhosh fc Lambl Il979allbh picture of 
neutron star accretion, a thin disc is threaded by the stel- 
lar magnetic field. The torque on the star is the sum of 
the material torque -/V mat carried by material falling onto 
the star and the torque Nb due to magnetic stresses. The 
material torque -/V mat always acts to spin up the star, 
while the magnetic torque Nb may either spin up or spin 
down the star. The Ghosh-Lamb picture has been extended 
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in se veral directions to acc ommodate time-dependent ef- 
fects. lLovelace et alJ il995t) proposed that the magneto- 
sphere could alte rnate between configu r ations with and 
without outflows. iLi fc Wickramasinghd il998T> suggested 
that the maximum radius where the stellar magnetic field 
threads the disc is constrained to small and large values, 
with intermediate values excluded by an unstable feedback 
mechanism as t he magneto s pheric configuration alters the 
disc resistivity. Torkelsson ( 1998) considered the interac- 
tion between a disc dynam o and the stellar magnetic field. 
Ivan Kerkwiik et alJ lll998l) considered the possibility that 
the i nner disc flips over, du e to warping, to form a retrograde 
disc. lYi fc Vishniad (1999) suggested that the accretion flow 
changes abruptly from a geometrically thin and cool Kep- 
lerian flow to a geometrically thick and hot sub-Keplerian 
flow, accompan i ed by a change from spin-up to spin-down. 
ISpruit fc Taaml dl993h (ST93 hereafter) showed that the disc 
can become unstable, with mass accreted in bursts, based 
on the premise that accretion is enhanced by the Kelvin- 
Helmholtz instability (KHI) but inhib ited when the magne- 
tosph ere rotates faster than the disc illlarionov fc Sunvaevl 
1975). They compared the bursts predicted by their model 
to the type II bursts observed from the Rapid Burster 
MXB 1730-355. 

In this paper, we generalize the accretion model of ST93 
to include torque feedback onto the star and diffusive pen- 
etration of the stellar magnetic field into the inner edge of 
the accretion disc, balanced by inward advection. On the 
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basis of preliminary numerical simulations, we also neglect 
enhancement of accretion by the KHI. We present the new 
model in Section In Section we show that the disc- 
magnetosphere system has two stable torque states, for cer- 
tain combinations of the magnetic Prandtl number Pm and 
the ratio £ of the corotation radius to the Alfven radius. If 
the star is initially spinning up, in the absence of extrane- 
ous perturbations, the spin-up equilibrium eventually van- 
ishes and the star subsequently spins down. In Section|l] we 
show that the viscous torque in the mixing layer, omitted in 
the original ST93 model, is significant; when it is included, 
the region of parameter space for which the system has two 
stable torque states is reduced. In Section^] we present the 
results of preliminary numerical simulations which indicate 
that the transport of material across the disc-magnetosphere 
mixing layer due to the KHI is suppressed when the mix- 
ing layer is thin. In Section HJ we discuss the implications of 
the generalized model for the torque reversals observed in 
accreting X-ray pulsars. 



2 DISC-MAGNETOSPHERE COUPLING 
2.1 Magnetospheric radius 

We use cylindrical polar coordinates (r, (j>, z) throughout and 
assume axisymmetry. The magnetospheric radius r m i of the 
system is defined to include all material that corotates with 
the star (ST93). [Such a radius may not exist if the magnetic 
field lines ar e swep t back all the way to the stellar surface; 
Ghos h et alJ 1^977)]. Outside r m i, there is a thin, annular 
mixing layer r m i < r < r m 2 within which the angular veloc- 
ity of infalling material is brought from the Keplerian value 
^K(rm2) to the stellar velocity f2* by magnetic torques. Once 
the material enters the mixing layer, it is transported on to 
the star on a time-scale much shorter than r m i/r m i. The 
rate at which mass flows inward through a surface at a fixed 
radius r in the disc is 



M(r,t) = -2-Krv r (r,t)S(r,i), 



(1) 



where E is the surface mass density, and v r < is the radial 
drift velocity. In the steady state, M is constant, indepen- 
dent of r (by continuity) , and equals the accretion rate onto 
the star. Out of the steady state, the accretion rate onto the 
star is 



M'(t) = M(r m2 ,t) + 2nr m2 r m 2T,(r m2 ,t), 



(2) 



where the second term in represents mass swept up as 
r m 2 changes. 

We place an upper bound on the steady state value of 
r m i by noting that the torque required to maintain corota- 
tion in the region r < r m i cannot exceed that provided by 
the magnetic shear stress, yielding (ST93) 



r-mi, 



MO, 



7v\S(r n 



(3) 



where S(r), the magnetic stress at radius r, at the d isc sur- 
face (z = h), is given by jShaoiro fc Teukolskvlll98^ 



\S(r)\ = (2*)- 1 \B,{r)B+(r)\ M=h . 



(4) 



The toroidal field, generated by differential rotation, is lim- 
ited to < B z , because as B$ builds up, the magnetic 



pressure inflates the field lines outward until they open up, 
destroying the m agnetic link between the star and the disc 
l)Alvlll98,4ll98i l. Hence, we have 



\S(r)\= V B 2 (r) (4tt) 



2 -&/a ^ 
rj/j, r (4-irj 



(•») 



where B(r) is the strength of the dipole field in the absence 
of a disc, /x is the stellar dipole moment, and r\ ~ 1 encapsu- 
lates our ignorance of the true (time dependent) magnetic 
field modified by the disc. 

The system does not necessarily adjust to r m i = 
J"mi,max. The value of r m i depends on how rapidly the mag- 
netic field (which induces corotation) penetrates the disc by 
Ohmic diffusion. Let ds = r m 2 — r m i be the depth to which 
the magnetic field penetrates the disc beyond r m i, and let 
Aj be the change in specific angular momentum across the 
mixing layer, given by 



Aj = r ml f2* - r m2 K (r- m 2). 



(6) 



Define d c to be the critical depth to which the field would 
need to penetrate in order to change the specific angular 
momentum of inflowing material by Aj, obtained by setting 
MAj = -27rr^ 1 d c S(r m i). One then finds 



d c = 



M(r m2 )n» 
27rS(r ml ) 



' ml / V / g J 



-3/2 



(7) 



where the corotation radius is given by r c = (GAf/f^) 1 / 3 , 
and our sign convention is such that a positive magnetic 
stress S subtracts angular momentum from the disc, as in 
ST93. Since the magnetic stress acts in the direction to bring 
material crossing the mixing layer into corotation, one has 
sgn(5) = sgn(r c - r m i). 

In equilibrium, one has ds = d c . If the equilibrium is 
perturbed by increasing M, d c increases, inflowing material 
cannot be brought into corotation as it travels from r m 2 to 
r m i, and r m i must decrease by an amount w d c — ds on the 
local radial drift time w ju r (r m 2)| /ds, giving 

ds — d c 



dr 



|«r-(r m2 )| . 



(8) 



Similarly, if the equilibrium is perturbed by reducing M, 
d c decreases, inflowing material is brought into corotation 
before travelling the full distance from r m 2 to r m i, and 
r m i must increase by an amount w ds — d c in a time 
~ K-(r m2 )| /d s , as in JSJ. 

2.2 Mixing layer thickness 

The mixing layer thickness ds is determined by a compe- 
tition between outward diffusion and inward advection. Let 
vb be the speed at which the magnetic field penetrates into 
the disc, as measured in a frame where v r = 0. One finds 
vb ~ va/dB, as for any diffusive process, where vb is the 
magnetic diffusivity of the disc material. In the observers 
frame, where the magnetospheric radius is changing and 
matter flows radially inward, we have vb = r m 2 + \ v r (rm2)\, 
and hence 



r m 2 



VB 



K(r m 2)| ■ 



(9) 



It might appear, from that we are erring by modelling 
a diffusive process as an advective one. This is not so, as 
is clear from the special case f m i = 0, Vr(r m 2) ~ 0, where 
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(|5Jl reduces to d,B = i^s/ds and hence ds = (2^st) 1//2 for 
= 0) = 0, as expected for diffusion. 

Comparing our model with ST93, we note the following. 
ST93 set r m i to the largest value for which the corotation of 
material within r m i can be maintained by magnetic stresses. 
Subsequently, they calculate the mixing layer thickness by 
requiring that disc material be brought into corotation by 
magnetic stresses before reaching r m i. In our model, the 
mixing layer thickness is set by the advection-diffusion bal- 
ance, and r m i grows or shrinks depending on whether or 
not material can be brought into corotation as it crosses the 
mixing layer, as described in Section l2.ll ST93 allow for r m i 
to decrease if the velocity shear across the mixing layer ex- 
ceeds the radial drift velocity, arguing that in this case the 
KHI transports material across the boundary layer faster 
than it can be brought into corotation. In our model, mate- 
rial is transported across the mixing layer at the radial drift 
velocity. This is motivated by numerical simulations of the 
KHI (Sectionjnj showing that mass transport is inhibited in 
a thin mixing layer. 

In both our model and ST93, magnetic stresses act only 
on the mixing layer, so the magnitude of the magnetic torque 
on the disc is smaller than in models of star-disc interac- 
tions where magnetic st resses act on the entire disc, such as 
IRappaport et all J2004D (RFS04 henceforth). For an order- 
of-magnitude comparison, one can take the representative 
values cLb = 0.01r c , r m i = l.lr c , and integrate our assumed 
form of the magnetic stress ISfV)! = /i 2 r _4 (47r) _1 over the 
mixing layer, to determine the approximate magnitude of 
the magnetic torque on the mixing layer in our model. One 
can similarly calculate the magnetic torque on the disc in the 
RFS04 model, by integrating the assumed form of the mag- 
netic stress in RFS04, \S{r)\ = /* 2 r~ 4 (2ir) _1 [l - fi K (r)/fi„], 
over the entire disc outside r c . We find that the magnetic 
torque on the mixing layer in our model is of order 0.03 
times the magnetic torque on the disc in RFS04, assuming 
identical stellar field strengths and spin frequencies in the 
two models. For an alternative perspective, one can take the 
assumed form of the magnetic stress in RFS04 and compute 
the ratio of the magnetic torque integrated over the entire 
disc outside r c to the torque integrated only over the region 
l.lr c < r < l.llr c . We find that these torques are in the 
approximate ratio 1:0.008. 



2.3 Surface density of the disc 

The surface density of the disc outside r m 2 evolves according 
to the thin disc equation 



as 

dt 



3d_ 

r dr 



1/2 



d_ 

dr 



(^Er 1/2 ) 



(10) 



where v is the kinematic viscosity. We neglect the possibility 
that the magn etic pressure in the mixing layer causes the 
disc to thicken dWandll987l) . The outer boundary condition 
is 



lim M(r, t) = M a 



(11) 



where Moo is some value set by the details of Roche lobe 
overflow. At the inner boundary, if r m 2 > r c , the mag- 
netosphere presents a centrifugal barrier and one expects 
S(r m 2) to be larger than if r m 2 < r c [ST93; see also 



RFS04], so we choose a boundary condition of the form 
j/E(r m 2) tx exp[a(r m 2/r c — 1)], with a > 0. We fix the con- 
stant of proportionality by considering steady state solu- 
tions, which have the form 



z/E : 



M 
3^ 



1/2 



(12) 



with /3 < 1 iFrank et al JI19S33) . If r m 2 = r c, the magneto- 
sphere corotates with the mixing layer, so the disc structure 
does not depend on the location of r m 2 (ST93) and 1121 
implies = and i/T, = M/(37r). We therefore choose the 
boundary condition at the inner edge to be 



Af(r m2 ) , , , 
fE(r m 2j = — expo(r m 2/rc — 1) • 

Comparison of 11211 and 11311 shows that 
/3 = 1 - exp[a(r m2 /r c - 1)] 



(13) 



(14) 



in steady state. We set a = 1 in this paper, except where 
specified otherwise. We show in Section[3]and Section^Jthat 
our results remain qualitatively unchanged for < a < 2. 
The boundary condition 11311 differs slightly from that in 
ST93, where j/E(r m 2) oc exp[a((r m 2/r c ) 3,/2 — 1)] and there is 
also a power law dependence of uT,(r m 2) on the mass accre- 
tion rate. 

The viscosity v in ST93 is constant in time and de- 
creases as a power law in radius. In our model, we assume 
that v is constant in time and independent of radius for sim- 
plicity. Later (Section 13.51 . we show that including a power 
law dependence of viscosity on radius leaves our results qual- 
itatively unchanged. 

2.4 Stellar spin frequency 

The stellar spin frequency, f2* , evolves according to the equa- 
tion of motion 



m, = N, 

where I is the star's moment of inertia and 
N = M'n K (r m2 )ri 2 + 2V„(r m2 ) 



(15) 



(16) 



is the torque on the star, set equal to the rate of change 
of angular momentum inside the surface r = r m 2. There 
is no magnetic term in I16H since, by assumption, there is 
no magnetic stress on the disc for r > r m 2. The first term 
in 1161 1 is the material torque and the second term is the 
viscous torque, given by 

(17) 



N v = 27rrz/Er 



dr 



In the steady state, 1101 reduces to 



(18) 



From 1141 and 1181 . the steady state torque on the star 
is positive (negative) when r m 2 is less than (greater than) 
r c . Away from the steady state, this is no longer necessar- 
ily true, since 1141 and 1181 hold true only in the steady 
state. [As matter accretes onto the star, the star's moment 
of inertia changes slightly, but this is a second order effect; 
iGhosh et alJ l)l977l) .] We show that long-term changes in £1* 
alter the equilibria of the star-disc system (Section I3.4|l : in 
contrast, ST93 held f2* constant. 
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2.5 Order-of-magnitude estimates of physical 
quantities 

In this section, we estimate the key physical parameters in 
the problem for typical systems. Let Mi. 4, .Z45, P10, and M-g 
denote the mass, moment of inertia, spin period and accre- 
tion rate of a neutron star, normalised by M = 1.4 M Q , 
I = 10 45 gem 2 , P spin = 10 s, and M = 10" 9 Moyr" 1 re- 
spectively. 

The corotation radius, defined by Q* = Qk^c), is given 

by 



d B = v B j |i>r(r m 2)| = d c 



(26) 



r c = 7.78 x W 8 M^P^ 3 cm. 



(19) 



The characteristic torque ./V c h ar at the corotation radius is 

JW = M(GMr c ) 1/2 (20) 

= 2.40 x 10 34 M_ 9 M 1 2/ 4 3 P 1 1 /3 gem 2 s" 2 . (21) 

The corresponding rate of change of the spin frequency is 

fL, char /(27r) = 38.2 x IQ- 13 I^M- 9 Ml'lPll 3 Hzs -1 . (22) 

Pulsars exhibiting torque reversals spin up (or down) 
at a compar able rate, with IflJ / (2rr) = 7-60 Hz s _1 
jNelson et al]ll997l) . In Section we show that our model 
predicts spin up/down rates of this order of magnitude. 

The kinematic viscosity in a di sc is conventionally 
param etrized by the a parameter IShakura fc Sunvaevl 
ll973j). From the standard unmagnetised a-disc solution 



1 ll992h . 



^4.0x 10 12 a 4/5 M% 10 M-l /4 (r/W 9 cm) cm 2 s" 



(23) 



taking r ~ r c and 0.001 < a < 1.0, we conclude that v lies 
in the range 10 10 to 10 13 cm 2 s _1 . It is convenient to define 
the viscous time-scale corresponding to the length-scale r c , 



r„ = 2r 2 /(3v) 

= 4.06 x 10 4 M 2/3 p 4/3 ' 
with fi3 = i^/(10 13 



P 

.4 ^10 



(24) 
(25) 



s _1 ). Using the above estimates for 
v and r c , we find r„ in the range 1 to 10 3 day; this is the 
time-scale for to rque transitions in our model (Section 13.31 . 
By comparison, Torkelsson (1998) found r„ to be typically 
less than 0.5 day 

If the magnetic diffusivity is due to the same turbu- 
lence that is responsible for the k inematic viscosity [ e.g. th e 
magneto- rotational instability of Balbus & Hawlev (1991)], 
one expects the magnetic Prandtl number to satisfy Pm = 
vb/v ~ 1. By contrast, in the literature dealing with the 
diffusion of magnetic fields through discs, Pm is frequently 
considered as a free parameter in the r ange 10 ~ 2 < Pm < 
10 2 jReves-Ruiz fc SteDinskill99dlLubow et al.lll994D . The 
work in this paper is valid for Pm <C 1, since for higher Pm 
the mixing layer is not thin (Section 



3 TORQUE BISTABILITY 

In this section, we show that the model of the disc- 
magnetosphere interaction presented in Section possesses 
two stable equilibria, with opposite signs of torque, under 
certain conditions. In searching for equilibria, we note that 
the time-scale fi./fi, ~ 10 13 s is longer than the time-scale 
t„ for r m i and <1b to change. Referring to (SJ and ©, we 
observe that the condition for equilibrium is 



First, we obtain dimensionless expressions for ds, d c 
and r m i jmax in equilibrium. From I26L 113L Q, and the 
definition of the Prandtl number Pm = v B /v, we find 

2 r m2 

3 r c 



ds/r c 



• Pm exp[a(r m2 /r c 



(27) 



A thin mixing layer (ds <§C r c ) therefore requires a small 
Prandtl number (Section 12. 51 . From Q and (JSJ, we find 

d c /r c = C /2 V2 (r m i/r c ) 6 sgn(r c - r ml ) 



r_m2 

r. 



-3/2 



- 1 



where 

j- -2/7 / 
5 = V r c /rA 

is a dimensionless parameter, and 

4 \ V7 



', 2GMM 2 



/' 



(28) 



(29) 



(30) 



is the characteristic A lfven radius for spherical accretion 
l|Elsner fc Lamblll977^ . From 0, © and JUJ, we find 



3/10^-7/10 



(31) 



Anticipating that r m i ~ r c , we initially restrict our investi- 
gation to the regime £ ^ 2~ 3 / 10 , r m i imax ^ r c , so that the 
constraint r m i < r m i jmax is satisfied. 

We now find equilibria by plotting ds/r c and d c /r c 
against r m2 /r c and solving 1261 graphically. Equations 126H . 
1271 1 and 12811 always admit one solution for r m 2 > r c , which 
we denote A (see Fig. . This is a spin-down equilibrium, 
since the torque is negative whenever r m 2 > r c (Section l2.4H . 

3.1 Multiple equilibria 

Under certain conditions, two other equilibria become avail- 
able to the system. The top panel of Fig. shows that for 
£ = 0.35, Pm = 0.01, the curves of ds/r c and d c /r c as func- 
tions of r m 2/r c have only a single intersection, so the disc- 
magnetosphere system has only a single, spin-down, equi- 
librium, denoted A. However, setting £ = 0.40, Pm = 0.01, 
as in the bottom panel of Fig. one finds that the curves 
of ds/r c and d c /r c have three intersections, and hence the 
disc-magnetosphere system has three equilibria, denoted A, 



B, C. These equilibria have r 



1.08r c , r 



(B) 



0.90r c , 



and = 0.56r c , so A is a spin-down equilibrium while B 
and C are spin-up equilibria. Combining 1181 and 1201 . one 
obtains 



f2* = P (r m 2/r c ) 1/2 f2» jChar 



(32) 



giving flY 1 ' = -0.082 n», char , fl^ = 0.089 fi», char , and 
ni ' = 0.27f2, iC h ar . In Section [3.21 we show that equilibria 
A and C are stable to perturbations in r m i and r m 2, while 
equilibrium B is unstable. The system is thus bistable, with 
the two stable equilibria having opposite signs of torque. 

A single equilibrium trifurcates into three equilibria 
when the curves of ds/r c and d c /r c in Fig. are tangential 
at some r m 2 < r c . Solving for this condition numerically, we 
obtain Fig. [5] which shows the combinations of Pm and £ 
for which multiple equilibria exist. For a = 1, the critical 
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0.002 



0.2 0.4 0.6 0.8 1 1.2 

Tm2 / Ic 



ing layer. At r m 2 ~ r c , there is no velocity shear across the 
mixing layer, so d c ~ 0. As r m 2 is displaced from r c , the 
velocity shear increases and hence d c increases. For small 
enough r m 2, d c decreases again because the magnetic stress 
is proportional to r~f , so it becomes 'easy' for the magnetic 
field to bring material crossing the mixing layer into coro- 
tation. The shape of d c as a function of r m i suggests that 
even if ds was determined differently than in Section 12.21 
(e.g. with a more sophisticated model of diffusion, or modi- 
fied by the action of the KHI), dB(r m 2) and d c (r m 2) might 
still have multiple intersections and bistability might still 
occur. 




0.2 0.4 0.6 0.8 1 1.2 

Fm2 / r c 



Figure 1. Equilibrium states for the disc-magnetosphere system. 
The top panel shows ds / r c (solid curve) and d c /r c (dashed curve) 
as functions of r m 2/r c , for Pm = 0.01, £ = 0.35. There is only one 
equilibrium, marked A, corresponding to spin down. The bottom 
panel shows the same information but for Pm = 0.01, £ = 0.4. 
There are two stable equilibria, marked A and C, which corre- 
spond to spin-down and spin-up respectively, and one unstable 
spin-down equilibrium, marked B. 



multiple equilibria 




a = 0.0 
a= 1.0 
a = 2.0 



single equilibrium 



0.01 0.02 0.03 0.04 0.05 0.06 
Pm 

Figure 2. Combinations of Prandtl number Pm and ratio 
£ of corotation radius to Alfven radius, for which the disc- 
magnetosphere system has a single equilibrium and multiple equi- 
libria. Curves for ^ a $J 2 are plotted. 



value £ cr it separating the single and multiple equilibria is 
well approximated by 

0.287 



Ccrit = 1.36Pm u 



(33) 



with a fractional error of < 1% in the range < Pm < 0.05. 

The existence of multiple equilibria can be understood 
as follows. One can think of d c as measuring how 'difficult' 
it is for magnetic stresses to change the angular velocity of 
infalling material from QK(r m 2) to fi, as it crosses the mix- 



3.2 Stability 

Suppose that an equilibrium is perturbed by slightly in- 
creasing r m i. From (|HJ, if ds — d c > 0, r m i will continue 
to increase, whereas if ds — d c < 0, r m i will tend back to 
the equilibrium value. This suggests that those equilibria for 
which ds — d c is a decreasing (increasing) function of r m 2 
are stable (unstable). Thus we expect equilibria A and C to 
be stable, and equilibrium B to be unstable. 

Elaborating on the above argument, we calculate the 
evolution of a perturbation about each equilibrium. As a 
(rough) first approximation to the equations of motion, we 
assume that changes in the surface density profile of the disc 
occur rapidly compared to changes in r m i, so that E(r, t) can 
be regarded as passing through a succession of equilibrium 
states given by 11121 . In reality, r m i and E change on simi- 
lar time-s cales; £ adjusts o ver a length-scale r m i in a time 
~ r mi/ v llFrank et alJll993 . p. 70), while the radial drift ve- 
locity in a thin disc 



-3 d 



/Er 1/2 ) 



(34) 



Er l/2 Q r 

together with JSJl, implies r m i/r m i ~ r^/i/. None the less, 
in our approximation, we find the following evolution equa- 
tions 



d /r m A 
T "&t V r c ) 

d /r m2 



r m 2 



exp[a(l — r m 2/r c )] 



(35) 



= ^Pm-^ — exp[a(l - r m2 /r c )] (36) 

3 ds r m 2 



is given 



by J2HJ, 



r^'], we write 



where the instantaneous value of d, 
given by (EH . and d_g = r m 2 — »" m i. 

Linearising about an equilibrium [r^^ 
the equations of motion in the form t„±; = AijXj, with 

x = [(r m i - r^f)/r c ,(rm2 - r^)/r c ] and A*, = dfi/dxj. 
The eigenvalues of Ay give the growth rate of perturbations. 
Obtaining analytic expressions for the growth rates is not 
tractable, so we calculate them numerically. For equilibria 
A, B and C of Fig. Q we find that the dominant modes 
of a perturbation have growth rates of -7.1, 3.8, and -4.7 
respectively, in units of r^ 1 , indicating that A and C are 
stable, while B is unstable. 

We emphasise that the above analysis provides an indi- 
cation of stability only, since it is based on the assumption 
that changes in the surface density profile of the disc are 
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Figure 3. Transitions between equilibria, for the system in the 
bottom panel of Fig. ^ If the system starts with r m 2 larger 
(smaller) than , it switches to equilibrium A (C). 

rapid compared to changes in the position of the magne- 
tosphere. To prove stability rigorously, one must to solve 
the disc PDEs subject to our chosen boundary conditions 
(e.g. ST93). However, this lies beyond the scope of this pa- 
per, especially in view of the uncertainties surrounding the 
boundary conditions themselves (e.g. the true, time depen- 
dent nature of the magnetic field in the vicinity of the mixing 
layer is unknown). 

3.3 Externally driven evolution 

How does the bistable system in Section loMl switch from one 
equilibrium to another? One possibility is that an extraneous 
perturbation, such as a change in Moo, pushes the system 
from A, through B (which is unstable), to C, in Fig. Q To 
explore this possibility, we integrate the equations of mo- 
tion 1351 1 and 136H . starting near the unstable equilibrium 
B. Typical results are shown in Fig. |3] One finds that the 
system tends to either of the stable equilibria, depending on 
the initial conditions. In each case, we set the initial value of 
cLb to the value corresponding to equilibrium B, but we set 
the initial value of r m 2 to be slightly larger (smaller) than 
r m2> whereupon the system tends to A (C). In each case, the 
transition between equilibria occurs on a time-scale ~ r v . 

3.4 Long-term evolution due to torque feedback 

The system can also switch from one equilibrium to another 
as fi* changes secularly in response to the disc torque. To 
explore this possibility, let us assume that the system is 
initially in a state of spin- up (equilibrium C). From I29II . 
£ oc r c oc Q^ 2 ^ 3 , so as the star spins up, £ decreases, the sys- 
tem moves from the upper to the lower region of parameter 
space in Fig. [5] and equilibria B and C vanish. The system 
then occupies the remaining equilibrium, A, and spins down. 
From (1311) . r m i, max /r c oc £ -7 / 10 oc OJJ 15 , so as the star spins 
down, r m i imax /r c decreases, until the state r m i = r m i, max is 
reached. Ultimately, r m i, max approaches r c , the torque tends 
to zero, and the system stops spinning down. 

Due to the dependences of £ and r m i jmax /r c on f2», fi» 
undergoes fractional changes of order unity in this evolution- 
ary scenario. For example, suppose Pm = 1, £ = 0.4, and 
initially Q* = 1 in arbitrary units. Then the star spins up 




-10 -5 5 10 

n 

Figure 4. Growth rates of perturbations about equilibria A, B 
and C, for the system in the bottom panel of Fig. Q as functions 
of the power law n governing the disc viscosity. As n varies, the 
growth rates change in magnitude but not in sign. 

to O* = 1.16, and subsequently spins down, with SI, — » 0.40 
as t — > oo. 

3.5 Radial viscosity gradient 

Thus far we have assumed that the viscosity is constant 
in time and independent of radius, whereas in ST93 v was 
allowed to vary as a power law in radius: 

v(r) = „(r c ) (r/r c ) n . (37) 

Adopting this prescription in our model, we find that the 
locations of the equilibria presented in Section 13.11 are un- 
changed, provided that we take Pm to be uniform (for sim- 
plicity). This is because the derivations of I27H and 128H do 
not rely upon the assumption of uniform v. With regard to 
the stability analysis in Section f3. 21 one finds that the evo- 
lution equations 13511 and I36H are each modified by a factor 
of (r m 2/r c )" on the right hand sides, and the v appearing in 
1241 1 is given by v{r c ). Numerical evaluation of the dominant 
linear growth rate (maximum eigenvalue of Aij), for a range 
of n, reveals that the growth rate changes in magnitude but 
not in sign (Figure [IJ, so the stability of each equilibrium is 
qualitatively unchanged from the case of uniform viscosity. 

More realistically, one might expect the viscosity v to 
vary as a function of the disc temperature and density, and 
therefore Pm may not be uniform. The detailed study of this 
scenario is beyond the scope of this paper. 

3.6 Wind torque 

For the low magnetic Prandtl numbers adopted in this study 
(Pm « 10 -2 ), significant bending of the magnetic field takes 
place, wh ich might lead to the formation of a centrifuga l 
wind (e.g. lLubow et al .11994 iRev es-Ruiz fc Stepinskll 996). 
accompanied by a magnetic torque that removes angular 
momentum from the disc. 

We can investigate the effect of the wind torque on our 
calculations with a toy model, in which we assume that a 
fraction 8 of the field lines threading the mixing layer are 
open and form a wind, and the magnetic stresses associated 
with these field lines remove angular momentum from the 
mixing layer (wind torque). The remaining fraction (1 — 5) of 
field lines threading the mixing layer are closed (connected 
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to the star) and can either add or subtract angular momen- 
tum from the mixing layer, as in Section [2.1 1 depending on 
the relative velocities of the star and disc. For simplicity, 
we assume that the toroidal field associated with the open 
field lines at the disc surface is equal in magnitude to the 
toroidal field associated with the closed field lines, though its 
direction may be different, giving ISwindl ~ ?7/x 2 r _6 (47r) _1 . 

In this crude picture, the wind torque modifies the equa- 
tion for d c , since d c is calculated by equating the magnetic 
torque on the mixing layer to the rate of change of angular 
momentum of material crossing the mixing layer. In the case 
r m i < r c , both the open and closed field lines subtract an- 
gular momentum from the mixing layer, and one finds that 
equation 12811 for d c is unchanged. In the case r m i > r c , the 
closed field lines add angular momentum to the mixing layer 
while the open field lines subtract angular momentum from 
the mixing layer, and one finds that 128H is modified by a 
factor of 1/(1 - 25) on the right hand side. For < 5 < 1/2, 
this has the effect of elevating the section of the curve of 
d c to the right of r c in Figure so that equilibrium A oc- 
curs closer to corotation. For 8 > 1/2, one has d c < for 
r m i > r c , so equilibrium A vanishes, and it is possible that 
there are no solutions for equilibria; the toy model breaks 
down. 

In our calculations, with or without the toy model for 
the wind, we have neglected magnetic torques on the disc 
outside the mixing layer. It is possible that there is a wind 
and an associated wind torque in this region. The wind 
torque might even dominate the viscous torque (see for ex- 
ample Pelletier 1992), and the omission of this torque is a 
limitation of our model. 



4 VISCOUS TORQUE AT THE MIXING 
LAYER 

The net torque exerted on the mixing layer induces the coro- 
tation of inflowing material and hence governs the evolution 
of the magnetosphere. It is the sum of magnetic and viscous 
components. The magnetic component results from the mag- 
netic stresses integrated over the top and bottom surfaces 
of the mixing layer (at \z\ = h), while the viscous compo- 
nent results from the mechanical shear stress integrated over 
the surface r = r m 2. In this section, we show that the vis- 
cous torque, neglected by ST93, is actually ~ r^/ds times 
larger than the magnetic torque. When the viscous torque is 
included in our model, only one stable equilibrium is avail- 
able to the disc-magnetosphere system under a wider range 
of conditions than otherwise, although this depends sensi- 
tively on the velocity profile near r m 2. 



4.1 Ratio of viscous and magnetic torques 

Material within r m i corotates with the star, so the viscous 
stress vanishes at the inner edge of the mixing layer. By 
assumption, the angular velocity profile of the disc is ap- 
proximately Keplerian beyond r m 2. The viscous torque N v 
acting on the mixing layer is therefore 

\dr)r= 
37Ti/E(r m2 ) f2 K (r m 2) r^ 2 - 



N„ 



27rr m2 ^2> n 



=>-m2 
2 



(38) 
(39) 



The negative sign indicates that the viscous torque sub- 
tracts angular momentum from material in the mixing layer. 
In the steady state, at r m2 , the surface density is given by 
^E(r m2 ) = M(l - /3)/(37r), implying 

\N v \=MSl K (r aa )rl 2 \l-P\. (40) 

Note that, for fixed M, |7V„| does not tend to zero in the 
limit v — > 0, since 14011 is independent of v. In this limit, 
from I34L v r oc v tends to zero, and £ oc M/v r oc v^ 1 
grows without bound, in such a way that N u oc i/E remains 
fixed. 

The magnetic torque Nb on the mixing layer is 



Nb = — 2nr m idBS(r m i)r m i. 



(41) 



From with r m i « r m i, max for order-of- magnitude pur- 



poses, we obtain MQ,, ~ 7rr m i 
\N B \ « 2r ml d B Km«. 
Dividing lHol by we find 

N v 



\S{r n 



implying 



(42) 



N B 



1 fl K (r m2 ) 
2 



r 2 

™m2 



n,. 



Tm2 



1-01: 



(43) 



2d B 

where we have used Q.K{r-m2) ~ fi* (close to corotation) and 
r m i w r m 2 (thin mixing layer). The expression (1431 shows 
that the viscous torque is much greater than the magnetic 
torque for ds <C J" m 2- Neglecting the viscous torque is not 
permissible, except for (3 close to unity. 



4.2 Equilibrium states 

In this section, we examine the modified equilibria of the 
disc-magnetosphere system after taking into account the vis- 
cous torque 13811 . Following a similar analysis to that pre- 
sented in Section r3.ll one finds that equation Ij27^ for ds/r c 
is still valid, but 128H needs to be modified as follows: 



d c /r c = £ r ' 2 V2 (r m i/r c ) 6 sgn(r c — r m i 
\r m i) \ r c J 



- 1 



(44) 



with B given by 114H . 

Examples of graphical solutions for equilibria, with the 
viscous torque included, are shown in Fig. |^| The top panel 
uses the same values of dimensionless parameters (Pm 
fi.f 1 1 . £ = 0.4) as the bottom panel of Fig. but clearly the 
shape of d c as a function of r m 2 is very different, and there is 
only one equilibrium. The bottom panel, with Pm = 0.001, 
£ = 0.4, shows a system with three equilibria; A and C 
are stable, with zero and spin-up torque respectively, while 
equilibrium B is unstable. In both panels, the discontinu- 
ity in d c arises because the magnetic torque changes sign 
at r m 2 = r c , while the viscous torque does not; in reality, 
the change in the magnetic torque occurs gradually over the 
domain r c — d_g < r m 2 < r c + ds- Fig. |H] shows the regions 
in Pm-£ parameter space for which there are multiple equi- 
libria. Note that, for a given £, the range of Pm giving rise 
to multiple equilibria is smaller compared to the case where 
the viscous torque is neglected (Fig. |5J . 

The viscous torque on the mixing layer depends on 
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Figure 5. Equilibrium states for the disc-magnetosphere system 
including the viscous torque at the mixing layer. The top panel 
shows ds/r c (solid curve) and d c /r c (dashed curve) as functions 
of r m 2/r c , for Pm = 0.01, § = 0.4. There is only one equilibrium, 
marked A, with zero net torque. The bottom panel shows the 
same information but for Pm = 0.001, £ = 0.4. There are two 
stable equilibria, marked A and C, corresponding to zero torque 
and spin-up respectively, and one unstable, spin-up equilibrium, 
marked B. 
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Figure 6. Regions in Pm-£ parameter space for which the sys- 
tem has a single equilibrium or multiple equilibria, with the vis- 
cous torque included. 

dQ./dr at r m 2 through 11381 1. The derivation of 1441 assumes 
that the disc velocity profile is Keplerian beyond r m 2, so that 
the viscous torque is given by 1391 . This assumption is ques- 
tionable, because viscous stresses may result in a nonkeple- 
rian velocity profile for a short distance beyond r m 2. Fig. [T] 
sketches two hypothetical velocity profiles in the vicinity of 
the mixing layer. In one case, the velocity profile is Kep- 
lerian outside r m 2, and the viscous torque on the mixing 
layer is negative, while in the other case the velocity profile 
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Figure 7. Schematic diagram of possible angular velocity profiles 
near the mixing layer (not to scale). Profile A is Keplerian be- 
yond r m 2 , and corresponds to a -ve torque on the boundary layer 
(dQ/dr < at r m 2). Profile B is non-Keplerian for a short dis- 
tance beyond r m 2, and corresponds to a +ve torque (d^l/dr > 
at r m 2). 



is non-Keplerian for a short distance beyond r m 2, and the 
viscous torque on the mixing layer is positive. Consequently, 
the system is very sensitive to the velocity profile near the 
mixing layer. ST93 did not address this issue, because the 
viscous torque on the mixing layer was neglected in their 
model. In the special case dQ./dr — at r — r m 2, the vis- 
cous torque on the mixing layer vanishes, and the analysis 
in Section |3] and ST93 is valid. 



5 TURBULENCE IN THE MIXING LAYER 

The interaction between the accretion disc and the stellar 
magnetic field at the disc-magnetosphere boundary is not 
well understood. The boundary is subject to magnetohy- 
drodynamic (MHD) instabilities, e.g. the magnetic inter- 
change instability (Mil) and the Kelvin-Helmholtz instabil- 
ity (KHI), which control the turbulent transport of material 
through the mixing layer ; into the ma gn etosphere, and onto 
the star iGhosh fc Lamblll97ct Il979allbt IScharlemannlll97gl : 
ISpruit fc Taamlll993l) . Such instabilities may also e xplain 
the quasi-periodic o scillations in X-ray sources feaanlll979l : 
iLi fc Naravanll2004ri . 

In this section, we assess the validity of the claim that 
the KHI transports material across the mixing layer with a 
speed (ST93) 



"KH = fml |^K(r m l) — fi» 



(45) 



To this end, we wrote a vortex-in-cell (VIC) code to per- 
form two dimensional (2D), inviscid, incompressible, hydro- 
dynamic simulations of the KHI (with B = and effective 
gravity g e g = 0) in an elongated domain, representing a thin 
mixing layer (ST93) . We emphasise that the simulations give 
some insight into the behaviour of the KHI but they do not 
faithfully simulate the global physics of the mixing layer. 
Ideally, in a full treatment, one would perform simulations 
in an unbounded domain and see how thick the mixing layer 
grows of its own accord, rather than enforcing a thin layer 
and testing what sort of flow is consistent with this assump- 
tion, as we do here. 
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5.1 Vortex-in-cell code 

The VIC algorithm <Liu fc Doorlvl2000l) discretises the vor- 
ticity field u) — V x v into point vortices, or vortons. The 
vortons are tracked as particles, while the velocity field 
v — V x A is obtained by solving Poisson's equation, 



0.05 



V^A = -u, 



(46) 



on a grid, where A is the vector potential for the velocity 
field. Vortons are advected by the velocity field; in turn, one 
can reconstruct the velocity field from the vorton distribu- 
tion. 

Our simulations are performed in a 2D box of dimen- 
sions < x < W = 1, -L/2 <y < L/2. The velocity field 
is tracked on a regular cartesian grid (x, y) tiled by n x x n y 
cells. The domain is periodic in one direction, 



;{x = 0,y, i) = v{x = l,y,t), 



(47) 



consistent with the cylindrical geometry of the annular mix- 
ing layer. In the other direction, representing the inner and 
outer edges of the annulus, we set 



v v (x,-L/2,t) = v y (x,L/2,t) = 0, 
w(x,-L/2,t) = w(x,L/2,t) =0. 



(48) 
(49) 



To be absolutely consistent with ST93, one would prefer 
to specify v x at \y\ = L/2. Doing so, however, requires 
the creation of new vortons at these boundaries, a feature 
that is hard to implement in our code. Initially, the fluids 
in the inner and outer halves of the annulus are approxi- 
mately counter-streaming, with v(x,y,0) — 0.5x for y < 
and v(x,y,0) = —0.5x for y > 0. Thus, one unit of time 
in our simulations is equal to the length of the simulation 
domain divided by the shear speed; physically, this is ap- 
proximately the spin period of the neutron star. To provide 
an initial perturbation from which the KHI can grow, each 
vorton is given a random initial y coordinate in the range 
-0.005 < y < 0.005, as in the top panel of Fig. |H1 



5.2 Transport speed 

The simulations track the root-mean-square and the maxi- 
mum of the y component of velocity, i) Sinlls and i> H , m ax, from 
which we can estimate the average and maximum rates at 
which material is transported across the mixing layer. We 
performed simulations for L = 1, 1/2, 1/4, 1/8, 1/16, 1/32, 
with n y = 64, n x = 64/L, and 5 x 10 3 vortons. Fig. [H]shows 
that the first vortical structures form on the smallest scales 
A, consistent with the linear KHI growth rate 7khi oc A -1 . 
These structures merge hierarchically until a single vortex 
occupies the simulation domain; energy cascades from small 
to large scales in 2D featchelorlll953l) . The top panel of 
Fig. [5] shows that i%, ma x rises initially, then decays within 
~ 20 crossing times, asymptoting to a value that depends 
upon W/L. The top panel of Fig. 1101 shows that higher 
aspect ratios W/L give a lower saturated value of i)j, mlx . 
The behaviour of iij, rms as a function of W/L is similar to 



%,max with Vy ^rms 



: 0.2vy t m&xL (bottom panels of Fig.|5|and 

Fig. EH . 

Our results imply that when the KHI is confined to 
an elongated domain, the maximum transport speed across 
the mixing layer is significantly less than the velocity shear. 
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Figure 8. KHI instability in a domain of dimensions W = 1.0, 
L = 1/8, at times t = 0, 0.125, 1, 10, 100 (top to bottom). Dia- 
monds represent vortons, and arrows show the velocity field. 



However, these results must be interpreted with caution. 
Our simulations are 2D; they neglect the Coriolis force, grav- 
ity, magnetic fields and compressibility; and the vorticity 
vanishes artificially at the boundaries. One would expect 
the inclusion of three-dimensionality, magnetic fields and 
compressibility to decrease i^.max furt her, since energ y cas- 
cades from large to small scales in 3D | |B^tcheIar|l953j) , com- 
press ibility tends to suppress the KHI ijWang fc Robertsonl 
1984) , and magnetic fields tend to disrupt large eddies 
jMalaeoli et alJll996Tl . In contrast, it is hard to predict the 
effect of different boundary conditions, the Coriolis force, 
and gravity (g c ff / 0). 



6 DISCUSSION AND CONCLUSIONS 

In this paper, we present a generalized model of the disc- 
magnetosphere interaction, which incorporates diffusion of 
the stellar magnetic field into the disc. We show that, for 
£ > Pm ' 3 , the system possesses two stable equilibria, cor- 
responding to spin-up and spin-down. (Our calculations in- 
dicate rather than prove stability, since they are based on 
the assumption that changes in the surface density profile 
of the disc are rapid compared to changes in the position 
of the magnetosphere). We suggest that transitions between 
stable equilibria can be induced, on the viscous time-scale 
r„ , by changes in the mass accretion rate at the outer edge of 
the disc (extraneous) or changes in the spin frequency of the 
star (torque feedback) . We also show that the viscous torque 
N v is generally much larger than the magnetic torque Nb; 
when it is included, the dynamics of the disc-magnetosphere 
system are very sensitive to the velocity profile at r m 2, but 
typically there are fewer combinations of Pm and £ for which 
the system is bistable. Preliminary numerical simulations of 
the KHI indicate that the transport speed across the mix- 
ing layer between the disc and magnetosphere is less than 
the shear speed when the layer is thin. In our model, mate- 
rial is transported across the mixing layer at the radial drift 
velocity. 

Given the presence of two stable equilibria in our model, 
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Figure 9. Maximum transport speed t) y , ma x (top panel) and 
root-mean-square transport speed v y>vnlE (bottom panel) due to 
the KHI across an elongated domain, versus time. One unit of 
time is approximately one stellar spin period. In each panel, 
curves are plotted for W/L = 1, 2, 4, 8, 16, 32 (top to bottom). 




Figure 10. The mean of t)j, m ai (top panel) and Bj^ms (bottom 
panel) over the last 5 time units of each simulation, as a function 
of W/L. One unit of time is approximately one stellar spin period. 



it is interesting to ask whether we can dra w any paral- 
lels w ith the phenomenon of torque reversals iNelson et alJ 
Il997ft . The X-ray pulsars GX 1+4, 4U-1626-67, Cen X-3 
and OAO 1657-415 have all been observed to alternate be- 
tween episodes of spin-up and spin-down, such that the time- 
scale for transitions between episodes (days) is much shorter 
than the the duration of individual torque episodes (weeks 
or years). A successful model for torque reversals should 
explain the observed properties: (i) Transitions never oc- 
cur between torques of the same sign (but different magni- 
tudes), (ii) The torque ranges from 0.2iVchar for 4U-1626-67 
to 2.3iV char for Cen X-3, assuming M = 10~ 9 M yr~\ (iii) 
The transition time-scale is much shorter than the time-scale 
for a sustained torque episode, (iv) The torque and luminos- 
ity are anti-correlated in the spin-down of GX 1+4. Evalu- 
ating our model against these criteria, we see the following: 
(i) The two stable equilibria have opposite signs of torque, 



so there are never transitions between torques of the same 
sign, (ii) The magnitude of the torque is ~ 0.1-/V c har- (hi) 
Transitions between equilibria occur over ~ 1 day (t„) for 
large values of viscosity {a. ~ 1), as observed, but are slow 
(~ 50 day) for a ~ 0.01 (iv) The anti-correlated torque and 
X-ray luminosity in GX 1+4 is not explained by our model. 
Also, our model predicts that, in the absence of extraneous 
perturbations, a star that is spinning up initially eventu- 
ally transitions to spin down, then evolves toward a state of 
zero torque, i.e. the model does not exhibit repeated torque 
reversals in its current, simple form. 

The model we have proposed, based on ST93, omits 
several effects which may be important. First, the disc 
may generate its own magnetic field via an MHD dynamo, 
driven by t he shear in the m ixing layer (as at the solar 
tachocline). iTorkelssonl dl99gft proposes dynamo action as 
a mechanism for torque reversals. Second, material may be 
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transported off the disc to form an outflow or be funnelled 
onto the stellar poles. For the low magnetic Prandtl num- 
bers adopted in our study (Pm w 10~ 2 ), significant bending 
of the field lin es takes place which might l e ad to a cen- 
trifugal wind jRev es-Ruiz fc Stepinskil fl99rj | Lubow et alJ 
1994). In the disc outside the mixing layer, the magnetic 
torque due to such a wind might dominate the viscou s 
torque felandford fe Pavnell982l:lPelletier fc PudritJl99l . 
whereas our model neglects magnetic torques outside of the 
mixing layer. Third, the mixing layer may not be thin in the 
radial direction. For example, even if the magnetic field does 
not penetrate far into the disc, the KHI may cause the mix- 
ing layer to thicken. Fourth, we assume the neutron star to 
be an aligned rotator, whereas the disc-magne tosphere sys - 
tem behaves differently for an oblique rotator JWandll997ft . 
Finally, our model does not solve for the disc structure in 
detail. Ideally, one would combine a detaile d solution for the 
disk structure (e.g. iRappaport et af]l2004h with a detailed 
model for the location of the inner disk radius (which in- 
cludes magnetic diffusion among other effects). 
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